*****************************************************************************
 *** Replication file for the paper:
 *** "You go, I stay: intrahousehold evacuation behavior upon a disaster"
 *** by Elias Maombi, Elie Lunanga, Nik Stoop, Marijke Verpoorten
 *** This version: March 2025
******************************************************************************* 
	
clear all
set more off
set scheme lean1
graph set window fontface "Latin Modern Roman"
	
********************************************************************************************
*** set working directory to where the database is saved
********************************************************************************************

cap cd "SET YOUR DIRECTORY"
	
********************************************************************************************
*** TABLES 
********************************************************************************************

*** Table 1: Overview of evacuation studies ==> based on literature review, see manuscript 


*** Table 2: Evacuation statuts of households and individuals 	
	use "data_volcano.dta", clear 
	
	* household evacuation status
	sum all_stay all_flee split  [iweight=ipw] if respondent==1 
	
	* individual evacuation status 
	sum fled [iweight=ipw] // all individuals
	sum fled [iweight=ipw] if respondent==1 // survey respondents
	sum fled [iweight=ipw] if hh_chief==1 // household heads
	
	* in partially evacuated households
	sum total_stay [iweight=ipw] if split==1 & respondent==1 
	sum prop_stay  [iweight=ipw] if split==1 & respondent==1 
	
	* household size
	sum hh_size  [iweight=ipw] if respondent==1


	
*** Table 3. Descriptive statistics of the explanatory variables 
	use "data_volcano.dta", clear 
	
	* Panel A 
	sum distance_lava probable_lava chief_age head_female tot_workingage total_chldren totage61_more total_malade proprietaire owns_business good_qual voitures family_propose [iweight=ipw] if respondent==1

	* Panel B  
	sum female child age61_more malade hh_chief [iweight=ipw] if ingoma==1
	
	

*** Table 4. Which households split? Bivariate analysis.

	* comparing "all members evacuated" to "some evacuated, some stayed" 
	use "data_volcano.dta", clear 
	keep if respondent==1 & all_stay==0
		
	reg log_distance split [pweight=ipw]
	reg probable_lava split [pweight=ipw]
	reg chief_age split [pweight=ipw]
	reg head_female split [pweight=ipw]
	reg tot_workingage split [pweight=ipw]
	reg total_chldren split [pweight=ipw]
	reg totage61_more split [pweight=ipw]
	reg total_malade split [pweight=ipw]
	reg proprietaire split [pweight=ipw]
	reg owns_business split [pweight=ipw]
	reg good_qual split [pweight=ipw]
	reg voitures split [pweight=ipw]
	reg family_propose split [pweight=ipw]
	
	
	* comparing "all members stayed" to "some evacuated, some stayed" 
	use "data_volcano.dta", clear
	keep if respondent==1 & all_flee==0
	
	reg log_distance split [pweight=ipw]
	reg probable_lava split [pweight=ipw]
	reg chief_age split [pweight=ipw]
	reg head_female split [pweight=ipw]
	reg tot_workingage split [pweight=ipw]
	reg total_chldren split [pweight=ipw]
	reg totage61_more split [pweight=ipw]
	reg total_malade split [pweight=ipw]
	reg proprietaire split [pweight=ipw]
	reg owns_business split [pweight=ipw]
	reg good_qual split [pweight=ipw]
	reg voitures split [pweight=ipw]
	reg family_propose split [pweight=ipw]


*** Table 5. Who evacuates in partially evacuated households? 
	use "data_volcano.dta", clear 
    keep if ingoma==1 & split==1

	eststo clear
	eststo: reg fled female child age61_more hh_chief malade i.evac_zone i.id [pweight=ipw], cluster(id) robust 
    eststo: reg fled female child age61_more hh_chief malade child_female i.evac_zone i.id [pweight=ipw], cluster(id) robust	
	eststo: reg fled female child age61_more hh_chief malade age61_female i.evac_zone i.id [pweight=ipw], cluster(id) robust	
	eststo: reg fled female child age61_more hh_chief malade chief_female i.evac_zone i.id [pweight=ipw], cluster(id) robust
	esttab using "table5.rtf", keep(female child age61_more hh_chief malade child_female age61_female chief_female) replace r2 nonotes nogaps



********************************************************************************************
*** Figures 
********************************************************************************************


*** Figure 1. Study area ==> maps, see manuscript 

*** Figure 2. Location of sample households and lava flow ==> map, see manuscript 

*** Figure 3. Determinants of household evacuation status, marginal effects 

	use "data_volcano.dta", clear
	keep if respondent==1
   	
	* Multinomial logit regression 
	// these estimated coefficients are presented in Table A5 in Appendix D
	local risk "log_distance probable_lava"
	local hh "chief_age head_female total_chldren tot_workingage totage61_more total_malade "
	local prop "proprietaire owns_business good_qual  voitures"
	local social "family_propose"
	
	mlogit hh_evac `risk' `hh' `prop' `social' evac_zone [pweight= ipw], baseoutcome(3) vce(robust) 
	estimates store m 
		
	* Marginal effects		
	// these marginal effects are presented in tabular form in Table A6 in Appendix D
	estimates restore m
	margins, dydx(*) predict (outcome(1)) post
	estimates store m1 
	
	estimates restore m
	margins, dydx(*) predict (outcome(2)) post
	estimates store m2 
	
	estimates restore m
	margins, dydx(*) predict (outcome(3)) post
	estimates store m3
	
	* Figure of marginal effects 
	// Figure 3 in the main manuscript 
	coefplot ///
	m1, bylabel({bf:Evacuation as a unit}) ///
	|| m3, bylabel({bf:Partial evacuation})  ///
	|| m2, bylabel ({bf:Staying as a unit}) ///
	xline(0) level(95) drop(_cons evac_zone) headings(log_distance="{bf:Risk}" chief_age="{bf:Household composition}" proprietaire="{bf:Property}" voitures="{bf:Car ownership}" family_propose="{bf:Social network}", labsize(small)) mlabel format(%9.2f) mlabposition(2) mlabgap(*3) mlabsize(vsmall)  msize(vsmall) m(d) ylabel(,labsize(small)) ciopts(lwidth(0.4)) byopts(compact rows(1)) subtitle(, size(small))	
	graph export "fig3.jpg", replace
	
	* Checking for multicollinearity 
	// These results are presented in Tables A7 and A8 in Appendix D
	reg hh_evac log_distance probable_lava chief_age head_female total_chldren tot_workingage totage61_more total_malade proprietaire owns_business good_qual  voitures  family_propose evac_zone [pweight= ipw], vce(robust) 
	vif
	
	collin log_distance probable_lava chief_age head_female total_chldren tot_workingage totage61_more total_malade proprietaire owns_business good_qual  voitures  family_propose evac_zone, corr
		   
	   
	
*** Figure 4. Number of individuals staying behind in partially evacuated households
	use "data_volcano.dta", clear
	keep if respondent==1 & split==1
 
	hist total_stay,percent xtitle(nr. household members staying) xlabel(1(1)12)
	graph export "fig4.jpg", replace
	 

     
*** Figure 5. Determinants of the nr. of evacuees in split households.
	use "data_volcano.dta", clear
	keep if respondent==1 & split==1

 	local risk "log_distance probable_lava "
	local hh "chief_age head_female tot_workingage total_chldren totage61_more total_malade"
	local prop "proprietaire owns_business good_qual voitures"
	local social "family_propose"
	
	* Tabular results: presented in Table A9 in Appendix E
	reg total_flee `risk' `hh' `prop' `social'  evac_zone [pweight=ipw], robust
	
	* Figure 5, presented in the main manuscript 
	coefplot, xline(0) level(95) drop(_cons evac_zone) headings(log_distance="{bf:Risk}" chief_age="{bf:Household composition}" proprietaire="{bf:Property}" voitures="{bf:Car ownership}" family_propose="{bf:Social network}") title("Nr. of members evacuating", size(small)) mlabel format(%9.2f) mlabposition(2) mlabgap(*1) mlabsize(vsmall) msize(vsmall) m(d) ylabel(,labsize(small)) ciopts(lwidth(0.5))	
	graph export "fig5.jpg", replace
	
	
	
*** Figure 6. Gender differences in evacuation rates within households that partially evacuated
	use "data_volcano.dta", clear
	keep if split==1
	graph hbar fled_female fled_male, over(age_group) ytitle(% of household members evacuating)  yvar(relabel(1"female" 2"male")) blabel(bar,format(%9.2f))
	graph export "fig6.jpg", replace
		
  
    
	
********************************************************************************************
*** APPENDIX Tables and Figures not yet created above
********************************************************************************************
	
*** Table A10: Reasons to stay
	
	use "data_volcano.dta", clear
	keep if all_flee==0 & respondent==1

	sum one_risk one_protect one_nowhere one_opport one_means one_sick one_young [iweight=ipw]
	sum one_risk one_protect one_nowhere one_opport one_means one_sick one_young [iweight=ipw] if split==0
	sum one_risk one_protect one_nowhere one_opport one_means one_sick one_young [iweight=ipw] if split==1
		

		
*** Table A11. Reasons to evacuate and means of transportation
	use "data_volcano.dta", clear
	keep if respondent==1
    
	sum reason_left_10 reason_left_11 reason_left_1 reason_left_6 reason_left_96 [iweight=ipw]
	tab final_transport [iweight=ipw]
	tab fuir_avant [iweight=ipw]		
	
	
	
*** Table A12. Gender, evacuation and risk aversion 
	use "data_volcano.dta", clear 
    keep if child==0 & ingoma==1

	reg fled female age61_more hh_chief malade log_distance probable_lava proprietaire owns_business good_qual voitures family_propose evac_zone [pweight=ipw]  if respondent==1, robust
	
	reg fled female age61_more hh_chief malade log_distance probable_lava proprietaire owns_business good_qual voitures family_propose evac_zone [pweight=ipw]  if respondent==1 & split==0, robust
	
	reg fled female age61_more hh_chief malade log_distance probable_lava proprietaire owns_business good_qual voitures family_propose  evac_zone [pweight=ipw]  if respondent==1 & split==1, robust
	
	reg fled female age61_more hh_chief malade log_distance probable_lava proprietaire owns_business good_qual voitures family_propose risk_aversion evac_zone [pweight=ipw]  if respondent==1 & split==1, robust


	
*** Table A13 in appendix. Gender and risk aversion
	use "data_volcano.dta", clear 
	keep if respondent==1
		
	reg risk_aversion female evac_zone [pweight=ipw], robust 	
	reg risk_aversion female hh_chief proprietaire owns_business good_qual log_distance evac_zone [pweight=ipw], robust
	
	

*** Table A14. Which household split ? Multinomial regression (eq 1) - without weights
	use "data_volcano.dta", clear 
	keep if respondent==1

	local risk "log_distance probable_lava"
	local hh "chief_age head_female tot_workingage total_chldren totage61_more total_malade"
	local prop "proprietaire owns_business good_qual voitures"
	local social "family_propose"

	 mlogit hh_evac `risk' `hh' `prop'  `social' evac_zone, baseoutcome(3) vce(robust) 
	

*** Table A15. Determinants of the number of evacuees in split HH - without weights
	use "data_volcano.dta", clear 
    keep if respondent==1 & split==1

	local risk "log_distance probable_lava "
	local hh "chief_age head_female tot_workingage total_chldren totage61_more total_malade"
	local prop "proprietaire owns_business good_qual voitures"
	local social "family_propose"

	reg total_flee `risk' `hh' `prop' `social'  evac_zone, robust

		
*** Table A16. Who, within split HH, evacuated? - without weights
	use "data_volcano.dta", clear 
    keep if ingoma==1 & split==1
	
 	local risk "log_distance probable_lava "
	local hh "chief_age head_female tot_workingage total_chldren totage61_more total_malade"
	local prop "proprietaire owns_business good_qual voitures"
	local social "family_propose"
		
	reg fled female child age61_more hh_chief malade evac_zone i.id, cluster(id) robust 
    reg fled female child age61_more hh_chief malade child_female evac_zone i.id, cluster(id) robust
	reg fled female child age61_more hh_chief malade age61_female  evac_zone i.id, cluster(id) robust	
	reg fled female child age61_more hh_chief malade chief_female evac_zone i.id, cluster(id) robust

 
	
*** Figure A2. Protecting property and quality of the construction. 
	
	use "data_volcano.dta", clear
	keep if respondent==1
	
		preserve
		* Some evacuated, some stayed
		keep if split==1
		
		egen group = group (constr_qual)
		gen protect_split=.
		sum group, meanonly 
		
		forvalues i= 1/`r(max)' {
		qui summ one_protect if `i'==group [iweight= ipw ]
		replace protect_split = r(mean) if `i'==group
		}
		drop group
		
		bysort constr_qual: gen x= _n
		keep if x==1
		keep constr_qual protect_split
		save protectsplit.dta, replace
		*
		restore
	  
		* All members stayed
		keep if all_stay==1
		
		egen group = group (constr_qual)
		gen protect_full=.
		sum group, meanonly 
		
		forvalues i= 1/`r(max)' {
		qui summ one_protect if `i'==group [iweight= ipw ]
		replace protect_full = r(mean) if `i'==group
		}
		drop group
		
		bysort constr_qual: gen x= _n
		keep if x==1
		keep constr_qual protect_full
	  
		merge 1:1 constr_qual using protectsplit.dta
		drop _merge
  
		label def qual 1 "very bad" 2 "bad" 3 "good" 4 "very good"
		label val constr_qual qual 
  
		twoway (line protect_full constr_qual, sort) (line protect_split constr_qual, sort), ylabel(0.4(0.1)0.8) xtitle("Quality of construction") ytitle("% of HH staying to protect property")  xlabel(1(1)4, angle(stdarrow) valuelabel) legend(order(1 "All members stayed" 2 "Some evacuated, some stayed"))
	


*** END ***	


	

	
	

	

	
	
	
	
	
	
	
	
	
	
	
	
